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ABSTRACT 


Earlier,  a novel  mathematical  model  of  buoyant  convection  in  an  enclosure  was 
developed.  The  nonlinear  equations  constituting  this  model  have  recently  been 
solved  by  finite  difference  methods  in  two  dimensions. 

In  this  paper  two  solutions,  obtained  in  special  cases,  to  the  model  equations 
are  presented.  For  both  cases  the  solutions  to  the  partial  differential  equa- 
tions and  to  the  finite  difference  equations  used  to  approximate  the  differen- 
tial equations  are  obtained  by  combinations  of  analytical  and  numerical  tech- 
niques. Agreement  between  the  exact  solutions  to  the  difference  equations 
described  in  this  paper  and  the  independently  obtained  numerical  solutions 
was  found  nearly  to  the  accuracy  specified  (usually  10“6)  for  an  iterative 
procedure  used  in  the  computational  scheme. 

The  first  solution  is  for  a time-dependent  irrotational , incompressible  flow  in 
an  enclosure  driven  by  sources  and  sinks  of  fluid  as  specified  by  the  heat 
source.  This  problem  arises  from  the  full  nonlinear  equations  with  boundary 
conditions,  in  continuous  or  discrete  form,  by  requiring  that  the  velocity 
field  be  irrotational  and  the  density  remain  constant. 

The  second  set  of  solutions  arises  when  several  other  simplifications  are  made 
to  the  equations.  The  density  is  taken  to  be  constant,  the  heating  is  assumed 
to  be  zero,  the  velocity  field  is  taken  to  be  two  dimensional  and  derivable 
from  a stream  function  only,  the  vorticity  is  taken  to  be  a constant,  and  the 
flow  is  independent  of  time.  These  solutions  are  used  to  determine  the  accu- 
racy with  which  the  code  described  in  Reference  2 solves  the  nonlinear  finite 
difference  equations  in  special  cases. 
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1.  Introduction 


Over  the  past  few  years,  the  National  Bureau  of  Standards  has  been  engaged  in  a 
research  project  to  develop,  starting  from  basic  conservation  laws,  a mathe- 
matical model  of  fire  development  within  a room.  Large  scale  convection  is  an 
essential  component  of  such  a model  because  this  fluid  motion  governs  the  smoke 
and  hot  gas  transport  within  a room  and  also  supplies  fresh  oxygen  to  the  fuel 
to  sustain  combustion.  Therefore,  development  of  a mathematical  model  of 
buoyant  convection  was  begun  as  a first  step  toward  a more  complete  room-fire 
model,  which  would  include  effects  of  combustion  chemistry,  radiation  and  smoke 
dynamics.  The  mathematical  model  for  convection,  the  partial  differential 
equations  and  boundary  conditions,  are  derived  in  Reference  1. 

Because  fluids  admit  a rich  variety  of  phenomena  and  because  the  equations  of 
fluid  dynamics  are  nonlinear,  it  is  difficult  to  obtain  solutions  to  these 
equations  except  in  very  special  cases.  It  is  the  nonlinear  nature  of  the 
equations  of  fluid  dynamics  that  makes  them  difficult  to  solve:  analytical 

tools  for  the  solution  of  nonlinear  problems  are  very  limited.  For  this  reason 
numerical  techniques  have  become  very  widely  used;  numerical  solution  of  non- 
linear equations  is  the  only  systematic  method  of  solving  them  under  a wide 
variety  of  conditions.  Yet,  if  these  numerical  solutions  are  to  be  trusted,  it 
is  essential  that  they  be  carefully  checked  to  determine  their  accuracy. 

It  is  the  purpose  of  this  paper  to  present  two  solutions,  obtained  in  special 
cases,  to  the  equations  derived  in  References  1 and  2.  Solutions  to  both  the 
partial  differential  equations  and  to  the  finite  difference  equations  are 
obtained  by  a combination  of  analytical  and  numerical  techniques  entirely  inde- 
pendent of  the  purely  numerical  procedure  described  in  Reference  2.  Agreement 
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between  the  exact  solution  to  the  difference  equations  described  in  this  paper 
and  the  numerical  solutions  described  in  Reference  2 was  found  nearly  to  the 
accuracy  specified  for  the  iterative  procedure  used  to  solve  the  pressure  equa- 
tion described  in  Reference  3 (usually  10"6).  On  the  other  hand,  comparison 
of  the  solution  to  the  difference  equations  and  of  the  corresponding  solution 
to  the  partial  differential  equation  determines  the  magnitude  of  the  discre- 
tization error  made  by  replacing  the  partial  differential  equations  by  finite 
difference  equations.  For  the  number  of  grid  points  used  in  most  of  the  cal- 
culations performed  to  date,  the  discretization  error  was  a few  percent.  In 
Reference  4 an  additional  analytical  solution  is  presented,  which  was  also 
used  to  assess  stability  and  accuracy  of  the  numerical  scheme. 

In  Section  2 the  partial  differential  equations  are  presented.  These  equations 
already  are  restricted  from  the  more  general  ones  presented  in  References  1 and 
2.  The  grid  is  also  shown  upon  which  the  discretized  variables  are  defined. 

The  first  solution,  presented  in  Section  3,  is  for  a fully  three-dimensional 
time-dependent  irrotational , incompressible  flow  in  an  enclosure  driven  by 
sources  and  sinks  of  fluid  as  specified  by  the  heat  source.  This  problem 
arises  from  the  full  nonlinear  equations  with  boundary  conditions,  in  con- 
tinuous or  discrete  form,  by  requiring  that  the  velocity  field  be  irrotational 
and  the  density  remain  constant.  The  resulting  nonlinear  partial  differential 
equations  and  difference  equations  can  then  be  solved.  The  solution  to  the 
difference  equations  is  used  to  determine,  in  detail,  the  accuracy  with  which 
the  code  described  in  Reference  2 solves  the  nonlinear  finite  difference  equa- 
tions in  this  special  case.  The  solution  to  the  corresponding  continuous 
problem  can  be  used  to  assess  the  discretization  errors  involved  as  a function 
of  mesh  size. 
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The  second  set  of  solutions  is  presented  and  discussed  in  Section  4.  This 
problem  arises  when  the  density  is  taken  to  be  constant,  the  heating  is  assumed 
to  be  zero,  the  velocity  field  is  taken  to  be  two  dimensional  and  derivable 
from  a stream  function  only,  the  vorticity  is  taken  to  be  a constant,  and  the 
flow  is  independent  of  time.  The  discretized  problem  can  be  solved  exactly, 
using  a combination  of  analytical  and  numerical  techniques  and  this  solution  has 
been  used  to  check  the  code  described  in  Reference  2. 

2.  Equations 

The  general  problem  considered  is  that  of  the  fluid  motion  in  a rectangular  en- 
closure produced  by  a specified  heat  source.  The  equations,  derived  in 

Reference  1 and  rewritten  in  more  convenient  form  for  numerical  integration  in 
Reference  2,  are  for  an  inviscid,  non-heat-conducting  perfect  gas.  In  the 
solutions  presented  in  this  paper,  density  variations  are  taken  to  be  small, 
and  buoyancy  effects  are  neglected.  Then  the  effects  of  density  variations  do 
not  alter  the  flow  to  a first  approximation. 

The  equations  describing  such  a flow  fluid  are  presented  below.  They  are  a 
simplified  version  of  Equations  (11)  of  Reference  2. 
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Q(x,y,z,t)  = QqX 


,p){T)  tanh  At  exp  [-$x(x-xc)‘ 


- 3y(y-yc)  “ Xz^ 


(2d) 


(2e) 


Po  is  a constant,  reference  density,  Q(x,y,z,t)  is  a heat  source  of  form  speci- 
fied in  (2e)  with  constants  3x»3y  and  x used  to  determine  its  spatial  extent  and 
A its  temporal  scale,  Qq  the  heat  source  magnitude,  and  P0(t)  is  a mean  pressure 
in  the  enclosure. 
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The  boundary  conditions  specify  that  there  will  be  no  outflow  or  inflow  at 
boundaries:  the  normal  velocity  is  zero.  Hence 

u = 0 or  v = 0 at  vertical  boundaries 

(3a) 

w = 0 at  horizontal  boundaries 

and  consequently  the  normal  derivative  of  pressure  at  boundaries  is  zero 

i-P  . 0 (3b) 

3n 

The  equations  can  be  made  dimensionless  by  introducing  a length  scale  equal  to  the 
height  of  the  enclosure,  a velocity  scale  determined  by  the  magnitude  of  the  diver- 
gence (or  the  heat  source)  specified  in  Equation  (2d),  the  ambient  density  pQ  and 
the  ambient  or  initial  pressure,  Poo.  When  this  is  done,  the  enclosure  becomes  one 
unit  high,  and  in  the  equations  above,  p , p*  and  Qq  can  be  taken  as  unity. 

In  Figure  1,  a two-dimensional  rectangular  enclosure  in  dimensionless  variables  is 
is  shown  together  with  a schematic  two-dimensional  representation  of  the  spatial 
grids  used  for  the  finite  difference  scheme.  The  grid  formed  from  solid  lines  re- 
presents the  basic  mesh  into  which  the  enclosure  is  divided:  in  general  there  are 

I mesh  cells  in  the  x-direction,  J mesh  cells  in  the  y-direction,  and  K mesh  cells 
in  the  z-direction. 

Upon  this  basic  mesh,  the  components  of  the  vector  (u,v,w)  and  the  component  of 
the  vector  vorticity  are  defined:  the  location  of  the  velocity  components  for  a 

sample  mesh  cell  are  also  shown  in  Figure  1. 
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A second  grid,  formed  by  joining  the  center  points  of  the  basic  grid  cells,  is 
that  upon  which  scalar  quantities  such  as  pressure  p are  defined.  In  Figure  1 
the  pressure  at  the  center  of  the  sample  mesh  cell  is  shown. 

The  following  discretely  evaluated  functions  will  denote  approximations  to  the 
corresponding  solutions  to  Equations  (9)  and  (10): 


u*1  = u(i5x,  (j-l/2)sy,  (k-l/2)sz,  nst) 

i jk 


v"  = v((i-l/2)sx,  jSy,  (k-l/2)sz,  nst) 
i jk 


w"  = w((i-l/2)6x,  (j-l/2)6y,  k6z,  nst) 
i jk 

(4) 


pn  s p( ( i -1/2 )6x,  ( j-1/2 )6y , (k-l/2)sz,  nst) 

i jk 


D°  s D ( ( i - 1/2 ) 6 x , (j-l/2)5y,  (k-l/2)«z,  n6t)  , 
i jk 


<j>n  = <j> ( (i-l/2)6x,  (j-l/2)5y,  (k-1/2 )6 z , nfit) 
l jk 


where  6x  = 1/(I-AR),  6y  = 1/ ( J • BR ) and  Sz  = 1/K  are  the  mesh  cell  sizes  in  the 
x- , y-  and  z-directions  respectively  and  where  St  is  the  time-step  size.  Such  a 
staggered  grid  is  commonly  used  for  multidimensional  finite  difference  integra- 
tions [5]. 
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In  the  following  two  sections,  solutions  to  nonlinear  equations  (1)  - (3)  under 
additional  restrictions  are  given.  Then  the  discretized  approximation  to  each 
restricted  set  of  partial  differential  equations  is  also  presented  and  solved. 

It  should  be  noted  that  each  of  the  discretizations  chosen  to  approximate  the  re- 
stricted set  of  partial  differential  equations  is  compatible  with  the  discretiza- 
tions chosen  in  Reference  2 for  the  more  general  partial  differential  equations 
describing  buoyant  convection  due  to  heating.  Because  of  this  compatibility,  the 
solutions  given  in  Sections  3 and  4 have  been  found  to  be  very  useful  for  compari- 
son with  the  numerical  results  obtained  from  the  general  algorithm  described  in 
Reference  2. 

3.  Irrotational  Flow  Driven  by  a Divergence 

The  nonlinear  problem  is  one  in  which  the  density  is  taken  to  be  constant  and  the 
velocity  is  assumed  to  be  derivable  from  a velocity  potential.  In  this  section, 
first  the  solution  to  the  partial  differential  equations  under  these  assumptions 
will  be  presented,  and  then  the  corresponding  solution  for  the  difference  equations 
is  given. 

A.  Continuous  Solution 

We  start  the  analysis  by  noting  generally  that  the  fluid  velocity  can  be  written 
as  the  a gradient  of  a scalar  potential,  plus  the  curl  of  a vector  potential 

+ + 

u = V4>  + v x ^ (5) 

For  an  irrotational  flow,  ^ = 0,  and  Equation  (Id)  becomes: 

V%  = D(x,y  ,z,t) 
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on  0 < x < 1/AR,  0 < y < 1/DR,  0 < z < 1.  The  boundary  conditions  are  found  from 

j\ 

Equation  (3)  to  be  — = 0 on  all  boundaries.  When  D(x,y,z,t)  is  separable,  the 


equation  for  <j>  is  separable,  and  <j>  can  be  determined  either  by  Fourier  series  or 
by  Green's  functions  and  complex-variable  techniques  in  the  two-dimensional  case. 
Then  the  spatial  distribution  of  <p  is  specified  by  the  spatial  distribution  of  D 
and  the  time  dependence  of  <p  is  equal  to  that  of  D.  The  velocity  field  is  found 
from 


The  pressure  field  is  found  from  the  dimensionless  version  of  Equation  (lb)  by 
noting  that  the  vorticity  is  zero.  With  the  density  constant  and  the  above  defi- 
nitions for  the  velocity  components  in  terms  of  the  potential.  Equation  (lb)  be- 
comes in  component  form 


3n 


u - 


(7) 


(8a) 


(8b) 


(8c) 


where  the  density  pQ  in  dimensionless  form  is  unity  and 


(9) 


Integration  of  Equations  (8)  yields 


(10) 
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The  function  g(t)  is  determined  from  the  condition 


1 1/AR  1/BR  , 

/ dz  / dx  / dyp(x,y,z,t)  = 0 (11) 

0 J0  o 

See  the  discussion  following  Equation  (8)  of  Reference  2 concerning  this  condi- 
tion. Integration  of  Equation  (10)  over  the  total  volume  0 < x < 1/AR, 

0 < y < 1/BR,  0 < z < 1 and  requiring  that 


/*  dz  /J/AR  dx  /^dy<t>(x,y,z,t)  = 0 (12) 

yields  for  g(t) 

1 1/AR  1/BR 
g(t ) = (AR-BR/2)  !q  dz  lQ  dx  fQ  dy 


The  pressure  field  is  then  found  from  the  potential  field  through  Equation  (10) 
where  g(t)  is  determined  from  Equation  (13).  Comparison  of  the  solution  obtained 
above  for  the  continuous  case  with  that  obtained  below  for  the  discretized  equa- 
tions then  provides  an  estimate  of  the  discretization  error  for  a specified  mesh 
size. 

B.  Discrete  Solution 

When  the  density  is  constant  and  the  discrete  velocity  components  are  derivable 
from  a discretized  potential,  a solution  to  the  difference  equations  can  be  ob- 
tained by  a procedure  analogous  to  that  used  above  to  obtain  a solution  to  the 
partial  differential  equations.  First,  reference  should  be  made  to  Figure  1 
where  the  discretized  potential  is  defined  at  the  center  of  a grid  cell.  Then, 
the  discretized  velocity  components  and  the  potential  are  related,  to  second 
order  accuracy,  by  the  relations 
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n 1 / n n 

u = — <p  - <(> 

ijk  5x  v i+l,j,k  ijk 


(14a) 


n 1 / n n 

v = — <t>  - 

ijk  Sy  Vi,j+l,k  ijk 


(14b) 


n 1 / n n 

w = - — ( <t>  -<() 

ijk  6z  \ ij,k+l  ijk 


(14c) 


Substitution  of  these  relations  into  the  discrete  version  of  Equation  (Id), 
yields  for  l<i<I,l<j<J,l<k<K, 


1/n  n n\l/n  n n 

~ — - ( d>  - 2d>  +d>  ] + — - d>  -2d>  + d> 

Sx2  \ i+l,jk  ijk  i-l,jk/  6y2  V i,j+l,k  ijk  i,j-l,k 


1/n  n n \ n 

+ — - / <p  - 2<p  + <p  | = D 

3_2  \ iJ»k+1  ijk  i j»k-l/  ijk 


(15) 


where  D is  the  discretized  version  Equations  (2d)  and  (2e), 
ijk 

The  discrete  boundary  conditions  are 


n n n n 

d>  = <t>  and  d.  = d>  for  OcjcJ  ,0<k<K 

0,jk  vl,jk  i+I , jk  vI,jk 


(16a) 


n n n n 

4>  ^ = <b  and  d>  = d>  for  0<i<I,0<k<K 

i,0, k i,l, k i,J+lk  yi,J,k 


(16b) 


n n n n 

and  <j> . = <j>.  for  0 < i < I , 0 < j < J 

ij,0  i j ,1  i J »K+1  l j ,K 


(16c) 
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n n n n n n 

where  d>  , d>  , <b  , <p  and  <j>.  . ^ , <j> 

0,jk  I+l,jk  i,0,k  i,J+l,k  ij,0  ij,K+l 

are  values  defined  at  the  center  of  fictitious  cells  adjacent  to  and  outside  the 

region  considered.  These  values  are  eliminated  from  the  problem  when  Equations  (16) 

are  substituted  into  Equation  (15).  Equation  (15)  and  boundary  conditions  (16)  can 

be  solved  using  fast  direct  methods  such  as  a three-dimensional  generalization  of  the 

n 

package  of  Swarztrauber  and  Sweet,  Reference  6.  When  D is  a separable  function  of 

ijk 

i,  j,  k and  n as  is  the  case  for  the  forms  chosen,  the  spatial  portion  of  <p  for  all 

ijk 

n can  be  determined  by  only  one  solution  of  the  linear  algebraic  system  represented  by 
Equations  (15)  and  (16).  The  discretized  velocity  field  is  then  determined  from 
Equations  (14). 

The  discretized  pressure  field  is  derived  from  the  discretized  version  of  Equations 
(lb).  Since  the  vorticity  is  zero  and  p0  = 1,  these  equations  can  be  written 


n+1  n 

u = U + 6 
ijk  ijk 
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26  x 
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kqi  jk  i 


2_ 


1 /nn  _ n ^ 
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(Ha) 
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v = v +6 

ijk  ijk  ( 2Sy 


( n 

Iqi  ,j+l,k. 


( n 
Vqijk 


1 / n n 

<$y  \Pi,j+l,k  Pijk 


(17b) 
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n+1 
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ijk 
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= W +6 
ijk 
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26z 


V,. 


2 — 


1 / n n ' 

6z  \Pij,k+l  Pijk 


(17c) 


where 


6 = 


St 


2st 


for  n = 1 
for  n > 1 


(18a) 
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for  n = 1 


i jk 


i jk 
n-1 

i 
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i jk 


n 

i 

i jk 
n-1 

f 

i jk 


w 


i jk 


w 


i jk 
n-1 

i 

i jk 


(18b) 


for  n > 1 


1 1 

and  where  q is  defined  by  Equation  (2a)  in  discretized  form, 
i jk 


~n 

and  U 

i jk 


1/2 


i-1 » jk , 


etc. 


The  upper  bracketed  values  in  Equations  (18)  are  the  ones  to  be  chosen  during  a first- 
order,  explicit  starting  time  step  whereas  the  lower  value  is  to  be  chosen  for  a leap- 
frog time  step.  In  this  form,  therefore.  Equations  (17)  and  (18)  can  be  used  to  re- 
present succinctly  both  first-order  and  leap-frog  time  steps.  It  is  important  to  note 
that  the  time  step  may  be  reduced  during  the  computation  to  avoid  violating  the 
stability  criterion. 


n 

6t  _<  max 

1<A<1  l<j<J 

1 <k<K 


+ 


'liiljkl  IVljlcI  liiljkl 

+ + 

- 6x  Sy  Sz 


-1 


(19) 


When  this  is  done,  the  time-marching  algorithm  is  restarted  with  the  values  of  the 
dependent  variables  at  the  last  successful  time  step  used  as  initial  conditions.  A 
first-order  time  step  is  then  taken  using  the  new  time  step  size  to  obtain  values  of 
the  dependent  variables  at  two  levels,  and  then  the  leap-frog  is  resumed. 
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We  substitute  definitions  (14)  into  Equations  (17)  and  define  a potential 


i jk 


n 

<t> . .. 

1 jk 


n-1 

♦ . .. 
ljk 


(20) 


as  in  Equation  (18)  to  take  account  of  the  different  form  of  Equation  (17) 
during  a first-order,  starting  time  step.  Then  Equations  (17)  become 


1 

6x 


n+1  n / n 


♦ i+ljk  ' $i+l,jk  \ i+1 , jk/  n 

+ i—  + n 

6 2 i+1 , jk 


r n+1  n / n \ 

<Mjk  " ^ i jk  j v i j k / t _n 
6 2 Pi jk 


= 0 


(21a) 


1 

«y 


n+1  n / n 

♦ 1 »j+l»k  ' $i  ,j+l  ,k  \qi,j+l,k/  n 

— + + p 

6 2 i,j+l,k 


r n+1  n / n 

jk  " $ijk  \Hijk)  n 

6 2 i jk 


= 0 


1 

6Z 


n+1  n / n \ 2 

^ i j » k+1  ’ $ij,k+l  VQij.k+l)  n 

+ A-.- /.  + pn-  ,• 

* 2 MlJ 


Pi  j,  k+1  ‘ 


’ n+1  n 
<Mjk  " $ijk 


n 

^ijk/  n 
+ Pi  jk 


= 0 


(21b) 


(21c) 
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The  difference  equations  are  satisfied  by  a function  independent  of  i,  j and  k 


(22) 


Equation  (22)  is  the  discretized  version  of  Equation  (10).  To  determine  g11,  we 
must  apply  the  discrete  analogue  of  Equation  (11) 


I J K n 

«x  «y  62  l l l p = 0 

i=l  j=l  k=l  'Jk 


and  require 


I J K n 

I l I = 0. 

i=l  j=l  k=l 


Then 


1 I J K IM-j-jk 

I6x  J6y  Kfiz  gn  = gn  = j J \ 6x  6y  6z  v 

AR-BR  i=l  j=l  k=l  ^ 


or 


n AR-BR  l J K 

gn  = __  J l l 6x  6y  6z 

2 i=l  j=l  k=l 


^i+1, jk  ~ *i-l,jkN 
V 2sx  / 


2 2 
*1»j+1.k  ^i>j-l,k\  + /*i,jk+l  ~ ^i j ,k-l\ 


26y 


26  z 


(23) 


(24) 


(25) 


(26) 
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where  we  have  used  the  relation 


2 


2 


(■ 


q. 


i jk 


n 


) 


2 


♦l+l, jk  ~ ♦i-l.Jk 
26  x 


(27) 


2 


+ 


*i j ,k+l  ~ (Hj,k-l 
^ 26  z 


n 


n 


Therefore,  once  the  discretized  potential  has  been  determined  from  Equation 


The  solution  procedure  outlined  above  has  been  carried  out,  namely  Equation 
(15)  was  solved  numerically  using  a three-dimensional  Poisson  solver  and  then  the 
velocity  field  was  determined  from  Eqs.  (14)  and  the  pressure  field  from  Equa- 
tion (22).  The  results  in  the  two  dimensional  case  were  compared  with  those 

computed  from  the  algorithm  described  in  Reference  2 when,  in  this  algorithm, 

~n  n 

the  perturbation  density  p and  the  vorticity  u were  specified  to  be  zero. 

i J i j 

The  parameter  in  the  pressure  solver  described  in  Reference  2,  which  determines 
the  accuracy  to  which  the  linear  algebraic  system  for  the  discretized  pressure  is 
solved,  was  taken  to  be  e = 10"6.  This  value  was  chosen  as  a reasonable  com- 
promise given  the  machine  accuracy  (36  bits)  of  the  NBS  UNIVAC  1100  and  the 
iterative  nature  of  the  pressure  solver.  Agreement  up  to  a few  tens  of  time 
steps,  between  dependent  variables  (both  components  of  velocity  and  pressure) 
determined  by  the  two  different  procedures,  was  found  to  a few  parts  in  10“ 6 . 
Accuracy  of  the  numerical  procedure  described  in  Reference  2 was  found  to 
degrade  slowly,  being  a few  parts  in  10“5  after  a few  hundred  time  steps. 

It  is  interesting  to  describe  the  physical  content  of  the  problem  discussed  in 
this  section.  When  the  density  is  taken  to  be  constant,  the  driving  function 
D(x,y,z,t)  can  no  longer  be  assumed  to  result  from  heating  and  cannot  be  inter- 


(15),  the  discretized  pressure  field  is  obtained  from  Equation  (22)  using 
Equation  (26)  to  determine  gn  and  using  Equation  (27)  for 
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preted  as  a source  of  specific  volume.  Rather,  D(x,y,z,t)  must  be  interpreted 
as  a distribution  of  sources  and  sinks  of  fluid  in  an  enclosure.  With  this 
interpretation  the  potential  flow,  with  no  outflow  or  inflow  at  the  boundaries 
of  the  enclosure,  has  some  interesting  features,  which  are  described  briefly 
bel ow. 

The  form  of  D(x,y,z,t)  used  for  these  calculations  is  defined  by  Equations  (2c) 
and  (2d)  with  Q(x,y,z,t)  given  by  Equation  (2e).  (For  the  solutions  shown  be- 
low, Q is  centered  along  floor  of  the  enclosure.)  Since  the  continuous  po- 
tential is  given  by  Equation  (6),  with  the  corresponding  velocities  given  by 
Equations  (7),  and,  since  Q and  therefore  D is  separable  in  time,  the  spatial 
dependence  of  the  potential  and  the  velocity  components  remains  the  same 
throughout  the  history  of  the  problem.  The  function  D(x,y,z,t)  represents 
sources  of  fluid  of  greatest  strength  at  the  center  of  the  "heat  source"  Q and 
decreasing  in  strength  as  the  distance  from  the  center  increases.  At  some 
distance  from  the  center,  the  sources  become  sinks  so  that  the  integral  of  the 
sources  (and  sinks)  over  the  enclosure  is  zero.  Therefore,  the  velocity  field 
is  one  describing  flow  from  the  "heat  source"  to  the  remainder  of  the  enclosure 
as  shown  in  Figure  2 for  all  times  during  the  history  of  the  flow. 

In  contrast,  the  pressure  changes  dramatically  during  the  flow.  During  the 
early  portion  of  the  flow,  all  nonlinear  manifestations  can  be  neglected,  and 
/ \ 8 <p 

in  Equation  (10),  p = - — . The  potential  <t>  will  behave  qualitatively  like 

3 1 

the  negative  of  the  source  function  D;  i.e.,  where  the  source  is  maximum  (at 
the  center  of  the  "heat  source")  the  potential  is  minimum,  and  the  pressure  is 
highest,  as  shown  in  Figure  3.  Hence  the  pressure  is  maximum  where  D is  and 
decreases  to  negative  values  a way  from  the  "heat  source"  Q. 
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Later  in  the  flow,  when  the  source  is  asymptoting  to  its  final  value,  the  non- 

linear  terms  dominate  and  — becomes  relatively  small  in  Equation  (10).  At 

3 1 

any  particular  time  g(t)  is  a constant,  and  Equation  (10)  becomes  approximately 


P = 


(28) 


The  term  in  square  brackets  is  the  magnitude  of  the  velocity  squared.  This 
equation  states  that  the  pressure  is  minimum  where  the  velocity  magnitude  is 
largest,  and  this  velocity  magnitude  is  largest  where  the  gradient  of  the  "heat 
source"  Q is  maximum.  Figure  4 shows  contours  of  constant  pressure  late  in  the 
flow  when  the  heat  source  is  asymptoting  to  its  final  value. 


4.  Two-Dimensional,  Steady-State  Flow  with  Vorticity 
A.  Continuous  Solution 


When  the  density  is  constant,  no  heat  source  is  present  and  we  restrict  the  flow 
to  be  two  dimensional,  then  the  continuity  equation,  Equation  (Id),  becomes 


— + - = 0 
3x  3y 

and  the  velocity  components 
ponent,  the  stream  function 

3^ 

u = — , v = - 

ay 

In  the  steady  state,  the  x- 
Equation  (lb),  become 


(29) 


are  derivable  from  a vector  potential  of  one  com- 
ip,  (see  Equation  (5)) 


dip 

3X 

and  y-  components  of  the  momentum  equation. 


(30) 
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(31) 


1 3 1 ap 

- — (uz  + V2)  - Vw  = - — — 

2 3X  p0  8 X 


1 d r ? 1 3P 

- — (u2  + V2)  + Uw  = - — — 

2 ay  p0  ay 


where  p0  is  constant  and  the  scalar  vorticity  is  defined  by 


av  a u .. 
OJ  = — - — = - V2]P 
ax  ay 


(32) 


Define 


H = p/Po  + 1/2  (u2  + v2) 


(33) 


and  substitute  for  u and  v from  Equation  (30)  into  Equations  (31).  Then 


aH  ai|> 

— + (o  — = 0 

ax  ax 


aH  dip 

— + o)  — = 0 

ay  ay 


If  we  now  assume  that  u is  a constant 


(34) 


0) 


= + b 


(35) 


then  Equations  (34)  become 


(H  + b\p ) = 0 
ax 


— ( H + bip ) = 0 

3y 
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or 

H + bi|i  = C = constant  (36) 

where  C is  evaluated,  for  example,  by  requiring  that  the  mean  value  of  p, 
integrated  over  the  enclosure  be  zero. 

Equation  (32)  becomes 

V2*  = - b (37) 

subject  to  the  impervious  wall  boundary  conditions,  which  become 


4>  = 0 at  x = 0,  1 
and  y = 0,  1 


(38) 


Therefore,  procedurally,  one  can  solve  the  Helmholtz  equation  (37)  subject  to 
boundary  conditions  (38)  to  obtain  the  stream  function  4>.  Then  Equations  (30) 
yield  the  velocities,  and  Equations  (33)  and  (36)  the  pressure. 


B.  Discrete  Solution 

As  in  the  continuous  case,  when  the  density  is  constant,  no  heat  source  is  pre- 
sent and  the  flow  is  two  dimensional,  the  continuity  equation  becomes 


un.  . un  , . vn.  _ vn  . . 

1J  1-1, J + 1J  1.J-1  _ 0 


fix  <5y 

and  the  discretized  velocity  components  are  obtained  from  the  stream  function 
as  follows 


(39) 


n 

u 

i j 


'P 


n 

i ,j-l 


5y 


n 


v 

ij 


n 

1-1  > J 


fix 


(40) 
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The  momentum  equations.  Equation  (31),  become  in  discretized  form 


1 1 
P0  <5x 


(41a) 


(41b) 


where 


n 

O) 

i j 


Sx 


sy 


(42a) 


ij 


1 

2 


(42b) 


(42c) 


Equations  (40)  imply  that 


n 1 / n 

w . = — r ( ^ 

6X2  \ 1+l.J 


n n \ 1/n  n n \ , % 

2*.  . + , . ) + - 2*..  + . . )(43) 

ij  i-l, J/  sy  \ !»J+1  iJ  i,J-l/ 
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and 


n 1 / n n 

u = — U - ^ 

u 25y  \i,j+l  i,J-l 

i j 


1/n  n 

26x  \i+l,j  ’N-IJ 


ij 


(44) 


Equations  (41)  can  be  written 


1 

6x 


1 

+ - 
2 


1 V 


1 / n \ 2 n 

P..,  ./ Pn  “ « q. . - P.  .Ip 


2 \ i+1 , j ) i+l» j 0 2 l ij 


ij 


nl/n  n \ n 1/n  n 

to  — • I \b  -ib  i+io  | ib  - ib 

ij  26 x V i+l,j  i-1 , j / i,j-l  26 x \ i+l,j-l  i-l,j-l 


(45a) 
= 0 


1 


1 

+ - 

2 


1/n  \ 2 n 1 

2 (qT,j+l  ) + Pi ,j+l/P°  ' 2 


- P. ./P, 
IJ  ' 


nl/n  n\n  1/n 

(0  I - lb  ) + 0)  (\h 

ij  25y  V i , j+1  i,j-l/  i-1, j 26y  V i-1, j+1 


n 

1-l.j-l 


(45b) 
= 0 


As  in  the  continuous  case,  we  take 


n 

to  = + b 

i j 


(46) 


and  note  that 


1 / n \ n , 

- q.  . + P.  7 p0 

2 \ ij  / ij  0 


(47) 


b / n n n n 

+ r I’P.  . + 4'.  , . + 4>  + ) = C = constant 

4 \ ij  i-l, J i,j-l  i-1 , j - 1 
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The  constant  in  this  equation  is  evaluated  by  requiring  that  the  mean  value  of 
n 

p over  all  mesh  points  be  zero  (to  make  it  compatible  with  the  general  de- 
i j 

fi nit  ion  of  this  quantity  in  the  algorithm  described  in  Reference  2). 

To  determine  a solution  to  the  nonlinear  difference  equations,  therefore,  the 
discretized  Helmholtz  equation  obtained  by  combining  Equations  (43)  and  (46) 


with  boundary  condi tons  that  i|».  . = 0 on  the  boundaries,  i = 0,  I and  j = 0,  J, 


is  solved  using  a software  package  such  as  the  one  developed  by  Swarztrauber 
and  Sweet  in  Reference  6.  From  this  solution  for  the  stream  function,  the 
velocities  are  determined  from  Equations  (40)  and  the  pressure  from  Equation 
(47).  This  solution  procedure  involves  analytical  and  numerical  techniques 
totally  independent  of  the  numerical  algorithm  described  in  Reference  2. 

The  steady-state  solution  to  the  difference  equations  was  used  to  test  the 
algorithm  of  Reference  2 by  using  the  values  as  initial  conditions  to  the  code 
(for  both  initial  time  levels)  with  the  heat  source  set  to  zero  and  the  per- 
turbation density  also  set  to  zero.  The  computer  code  was  then  allowed  to  take 
several  hundred  time  steps,  and  the  values  after  different  numbers  of  time 
steps  were  compared.  It  was  found  that  there  was  no  instability  in  the  algo- 
rithm. The  solution  replicated  itself  very  well,  with  a slow  degradation  in 
the  number  of  significant  figures  with  which  the  solution  was  able  to  duplicate 
itself  as  the  number  of  time  steps  increased.  Initially,  the  solution  after  a 
few  time  steps  agreed  with  the  exact  solution  to  a few  parts  in  the  sixth 
significant  figure,  as  expected  from  the  tolerance  set  on  the  pressure  solver. 


(49) 
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After  four  hundred  time  steps,  the  agreement  was  to  a few  parts  in  the  fourth 
significant  figure,  representing  a degradation  which  we  felt  was  reasonable. 

In  Figure  5 we  present  plots  of  the  horizontal  velocity  as  a function  of  ver- 
tical location  at  three  horizontal  positions  and  of  the  vertical  velocity  as  a 
function  of  horizontal  location  at  three  vertical  positions.  This  figure  shows 
a flow  field  which  is  rotating  in  a clockwise  direction  with  the  greatest 
velocities  at  the  edges  and  smaller  velocities  in  the  interior.  The  flow  field 
represents  "stirring  in  a rectangular  tea  cup"  with  no  dissipative  effects  of 
viscosity.  Horizontal  velocity  plots  at  a specified  horizontal  position  do  not 
agree  with  vertical  velocity  plots  at  the  corresponding  vertical  location  due 
to  details  in  the  graphics  package  used  to  display  the  data. 
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Figure  Captions 


Figure  1: 

A rectangular  enclosure  in  dimensionless  variables  and  the  a sche- 
matic representation  of  the  basic  mesh  into  which  the  enclosure 
is  divided. 

There  are  I cells  in  the  x-direction,  J cells  in  the  y-direction 
and  K cells  in  the  z-di recti  on.  Also  shown  is  a single  mesh  cell 
and  the  location  within  this  mesh  cell  of  the  discretely  defined 
dependent  variables,  velocity  components  u,  v and  w,  pressure  p 
and  velocity  potential  <f>. 

Figure  2: 

The  velocity  field  in  the  two-dimensional  case.  The  horizontal 
velocity  is  plotted  at  three  horizontal  locations  (in  each  case  as 
a function  of  the  vertical  coordinate).  The  vertical  velocity  is 
shown  at  three  vertical  locations  as  a function  in  each  case  of  the 
horizontal  coordinate.  The  flow  field  is  from  the  sources  to  the 
sinks;  the  spatial  pattern  does  not  change  with  time  for  this 
potential  flow. 

Figure  3: 

Contours  of  constant  pressure  early  in  the  calculation  of  the 
two-dimensional  potential  flow.  Contours  greater  than  the  mean 
pressure  are  shown  as  solid  lines,  while  those  below  the  mean  are 
dashed.  The  pressure  is  maximum  where  the  distribution  of  sources 
is  and  decreases  to  negative  values  away  from  the  "heat  source" 
(where  the  sinks  are). 

Figure  4: 

Contours  of  constant  pressure  late  in  the  calculation  of  the 
two-dimensional  potential  flow.  Contours  greater  than  the  mean 
pressure  are  shown  as  solid  lines,  while  those  below  mean  are 
dashed.  Note  that  the  pressure  is  negative  where  the  sources  are 
and  positive  away  from  the  sources  in  contrast  to  Figure  3.  The 
change  in  character  results  from  the  Bernoulli  effect  (see  the  text 
for  an  explanation). 

Figure  5: 

Horizontal  velocity  as  a function  of  vertical  location  at  three 
horizontal  positions  and  vertical  velocity  as  a function  of  hori- 
zontal location  at  three  vertical  positons.  Flow  field  is  rotating 
in  the  clockwise  direction  with  greatest  velocities  at  the  edges 
and  smaller  velocities  in  the  interior,  representing  a stirring  in 
the  rectangular  region. 
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